#******************************************************************************
#******************************************************************************

# import libraries

import math

# TODO: convert functions to use objects

from src import pysoleng as pse

import datetime as dt

#******************************************************************************
#******************************************************************************

def test_solar_position_four_corners():
    
    # https://gml.noaa.gov/grad/solcalc/
    
    list_tests = [
        # north america
        {'latitude':34.307143, # positive north
         'longitude':-94.394531, # negative west
         'local_time':dt.datetime.fromisoformat('2022-07-15 19:35:28-05:00'),
         'true_et_minutes':-6.05, # minutes
         'true_solar_declination':21.39,
         'true_solar_noon':dt.datetime.fromisoformat('2022-07-15 13:23:33-05:00'),
         'true_sunrise_time':dt.datetime.fromisoformat('2022-07-15 06:17:00-05:00'),
         'true_sunset_time':dt.datetime.fromisoformat('2022-07-15 20:30:00-05:00'),
         'true_solar_azimuth':289.45-180, # original: 289.45
         'true_solar_elevation':9.64},
        # eastern europe
        {'latitude':52.268157, # positive north
         'longitude':26.71875, # negative west
         'local_time':dt.datetime.fromisoformat('2022-11-01 15:49:34+03:00'),
         'true_et_minutes':16.47, # minutes
         'true_solar_declination':-14.52,
         'true_solar_noon':dt.datetime.fromisoformat('2022-11-01 12:56:41+03:00'),
         'true_sunrise_time':dt.datetime.fromisoformat('2022-11-01 08:08:00+03:00'),
         'true_sunset_time':dt.datetime.fromisoformat('2022-11-01 17:44:00+03:00'),
         'true_solar_azimuth':222.99-180, # original: 222.99
         'true_solar_elevation':13.56},
        # south america
        {'latitude':-44.33956, # positive north
         'longitude':-69.785156 , # negative west
         'local_time':dt.datetime.fromisoformat('2022-11-01 15:49:34-03:00'),
         'true_et_minutes':16.48, # minutes
         'true_solar_declination':-14.6,
         'true_solar_noon':dt.datetime.fromisoformat('2022-11-01 13:22:41-03:00'),
         'true_sunrise_time':dt.datetime.fromisoformat('2022-11-01 06:19:00-03:00'),
         'true_sunset_time':dt.datetime.fromisoformat('2022-11-01 20:27:00-03:00'),
         'true_solar_azimuth':302.02-180, # original: 302.02
         'true_solar_elevation':46.98},
        # east africa
        {'latitude':-28.30438, # positive north
         'longitude':30.8935546, # negative west
         'local_time':dt.datetime.fromisoformat('2022-11-01 15:49:34+02:00'),
         'true_et_minutes':16.47, # minutes
         'true_solar_declination':-14.53,
         'true_solar_noon':dt.datetime.fromisoformat('2022-11-01 11:39:59+02:00'),
         'true_sunrise_time':dt.datetime.fromisoformat('2022-11-01 05:04:00+02:00'),
         'true_sunset_time':dt.datetime.fromisoformat('2022-11-01 18:16:00+02:00'),
         'true_solar_azimuth':269.44-180, # original: 269.44
         'true_solar_elevation':30.95},
        ]
    
    for mydict in list_tests:
        
        # get location and local time objects
        
        mylocation = pse.timeplace.Location(
            coordinates=(mydict['latitude'], mydict['longitude']),
            crs_string='epsg:4326'
            )
        
        mytime = mydict['local_time']
        
        #**********************************************************************
    
        # declination
        
        abs_tol = 1.1
        
        day_year = pse.timeplace.day_of_the_year(mytime, False)
        
        dec_angle = pse.solar.declination_angle(day_year)
                        
        assert math.isclose(dec_angle,
                            mydict['true_solar_declination'],
                            abs_tol=abs_tol)
        
        dec_angle_spencer = pse.solar.declination_angle_spencer1971(day_year)
        
        abs_tol = 0.2
        
        assert math.isclose(dec_angle_spencer,
                            mydict['true_solar_declination'],
                            abs_tol=abs_tol)
        
        #**********************************************************************
    
        # equation time
        
        abs_tol = 0.4 # 0.05863998
        
        et = pse.solar.equation_of_time_minutes(day_year)
        
        assert math.isclose(et,
                            mydict['true_et_minutes'],
                            abs_tol=abs_tol)
        
        # equation time by spencer
        
        abs_tol = 0.2 # 0.05863998
        
        et_spencer = pse.solar.equation_of_time_minutes_spencer1971(day_year)
        
        assert math.isclose(et_spencer,
                            mydict['true_et_minutes'],
                            abs_tol=abs_tol)
        
        # equation time by smart
        
        abs_tol = 0.1 # 0.05863998
        
        JD = pse.solar.julian_day(mytime.year,
                                  mytime.month, 
                                  mytime.day)
        
        T = pse.solar.time_in_julian_centuries_since_j2000(JD)
        
        mean_solar_anomaly = pse.solar.mean_anomaly_sun(T)
        
        obliquity_ecliptic = pse.solar.mean_obliquity_ecliptic_IAU(T)
        
        mean_solar_longitude = pse.solar.mean_longitude_sun(T)
        
        earth_orbits_eccentricity = pse.solar.eccentricity_earths_orbit(T)
        
        et_smart = pse.solar.equation_of_time_radians_Smart(
            mean_solar_anomaly*math.pi/180, 
            earth_orbits_eccentricity, 
            mean_solar_longitude*math.pi/180, 
            obliquity_ecliptic*math.pi/180)
                
        assert math.isclose(et_smart*(180/math.pi)*60/15,
                            mydict['true_et_minutes'],
                            abs_tol=abs_tol)
        
        #**********************************************************************
                
        # apparent solar time
        
        ast = pse.solar.lst_to_ast(mytime, 
                                    mylocation, 
                                    #equation_time_minutes=et_spencer
                                    equation_time_minutes=(
                                        et_smart*180/math.pi)*60/15
                                    )
        
        solar_ha = pse.solar.hour_angle(ast)
        
        # solar altitude
        
        abs_tol = 1
        
        solar_alt = pse.solar.solar_altitude_angle(
            mylocation.latitude*math.pi/180, 
            dec_angle_spencer*math.pi/180, 
            solar_ha)
        
        assert math.isclose(solar_alt*180/math.pi,
                            mydict['true_solar_elevation'],
                            abs_tol=abs_tol)

        #**********************************************************************
        
        # zenith angle
        
        cos_zen_angle = pse.solar.cosine_zenith_angle(
            mylocation.latitude*math.pi/180,
            dec_angle_spencer*math.pi/180, 
            solar_ha)
        
        cos_inc_angle = pse.solar.cosine_incidence_angle(
            declination_angle=dec_angle_spencer*math.pi/180,
            latitude=mylocation.latitude*math.pi/180,
            plane_slope_horizontal=0, 
            plane_azimuth_south=0,
            hour_angle=solar_ha)
        
        abs_tol = 0.01
        
        assert math.isclose(cos_zen_angle,
                            cos_inc_angle,
                            abs_tol=abs_tol)
        
        #**********************************************************************
        
        # solar azimuth
        
        solar_az = pse.solar.solar_azimuth_angle(
            mylocation.latitude*math.pi/180, 
            dec_angle_spencer*math.pi/180, 
            math.acos(cos_zen_angle), 
            solar_ha)
        
        abs_tol = 1
        
        assert math.isclose(solar_az*180/math.pi,
                            mydict['true_solar_azimuth'],
                            abs_tol=abs_tol)
                        
        #**********************************************************************
                
        # sunrise hour angle
        
        time_tolerance = dt.timedelta(minutes=7)
            
        true_lst = mydict['true_sunrise_time']
        
        # hour angle -> ast 
        
        solar_ha_sunrise = pse.solar.sunrise_hour_angle(
            dec_angle_spencer*math.pi/180,  # dec_angle_spencer
            mylocation.latitude*math.pi/180)
        
        sunrise_ast = pse.solar.ast(solar_ha_sunrise,
                                    true_lst.year,
                                    true_lst.month,
                                    true_lst.day,
                                    tzinfo=true_lst.tzinfo)
        
        sunrise_lst = pse.solar.ast_to_lst(sunrise_ast, mylocation)
        
        # test upper bound
        
        assert sunrise_lst <= true_lst + time_tolerance
        
        assert sunrise_lst >= true_lst - time_tolerance
        
        # sunrise hour angle using meeus method
        
        time_tolerance = dt.timedelta(minutes=1)
        
        solar_ha_sunrise = pse.solar.sunrise_hour_angle_meeus(
            dec_angle_spencer*math.pi/180,  # dec_angle_spencer
            mylocation.latitude*math.pi/180)
        
        sunrise_ast = pse.solar.ast(solar_ha_sunrise,
                                    true_lst.year,
                                    true_lst.month,
                                    true_lst.day,
                                    tzinfo=true_lst.tzinfo)
        
        sunrise_lst = pse.solar.ast_to_lst(sunrise_ast, mylocation)
        
        # test upper bound
        
        assert sunrise_lst <= true_lst + time_tolerance
        
        assert sunrise_lst >= true_lst - time_tolerance
                
        #**********************************************************************
                
        # sunset hour angle
        
        time_tolerance = dt.timedelta(minutes=7)
            
        true_lst = mydict['true_sunset_time']
        
        # hour angle -> ast -> lst
        
        solar_ha_sunset = pse.solar.sunset_hour_angle(
            dec_angle_spencer*math.pi/180,  # dec_angle_spencer
            mylocation.latitude*math.pi/180)
        
        sunset_ast = pse.solar.ast(solar_ha_sunset,
                                   true_lst.year,
                                   true_lst.month,
                                   true_lst.day,
                                   tzinfo=true_lst.tzinfo)
        
        sunset_lst = pse.solar.ast_to_lst(sunset_ast, mylocation)
        
        # test upper bound
        
        assert sunset_lst <= true_lst + time_tolerance
        
        assert sunset_lst >= true_lst - time_tolerance
        
        # sunset hour angle using meeus method
        
        time_tolerance = dt.timedelta(minutes=2)
        
        solar_ha_sunset = pse.solar.sunset_hour_angle_meeus(
            dec_angle_spencer*math.pi/180,  # dec_angle_spencer
            mylocation.latitude*math.pi/180)
        
        sunset_ast = pse.solar.ast(solar_ha_sunset,
                                   true_lst.year,
                                   true_lst.month,
                                   true_lst.day,
                                   tzinfo=true_lst.tzinfo)
        
        sunset_lst = pse.solar.ast_to_lst(sunset_ast, mylocation)
        
        # test upper bound
        
        assert sunset_lst <= true_lst + time_tolerance
        
        assert sunset_lst >= true_lst - time_tolerance
        
        #**********************************************************************
        
        # solar noon
        
        # hour angle -> ast -> lst
        
        true_lst = mydict['true_solar_noon']
        
        time_tolerance = dt.timedelta(minutes=1)
        
        solar_ha_noon = 0
        
        solar_noon_ast = pse.solar.ast(solar_ha_noon, 
                                       year=true_lst.year,
                                       month=true_lst.month,
                                       day=true_lst.day,
                                       tzinfo=true_lst.tzinfo)
        
        solar_noon_lst = pse.solar.ast_to_lst(solar_noon_ast, mylocation)
        
        # test upper bound
        
        assert solar_noon_lst <= true_lst + time_tolerance
        
        assert solar_noon_lst >= true_lst - time_tolerance
        
    #**************************************************************************

#******************************************************************************
#******************************************************************************

def test_obliquity_ecliptic():
    
    # Source: Meeus, 1998; Example 22.a, page 148
    
    # Date: 1987 April 10th, 0^h TD
    
    #**************************************************************************
    
    JDE = pse.solar.julian_day(1987, 4, 10)
    
    JDE_true = 2446895.5
    
    abs_tol = 1e-6
    
    assert math.isclose(JDE, JDE_true, abs_tol=abs_tol)
        
    #**************************************************************************
    
    # time in julian centuries since j2000
    
    T = pse.solar.time_in_julian_centuries_since_j2000(JDE_true)
    
    T_true = -0.127296372348
    
    abs_tol = 1e-6
    
    assert math.isclose(T, T_true, abs_tol=abs_tol)
    
    #**************************************************************************
    
    # mean obliquity of the ecliptic (IAU)
    
    OE = pse.solar.mean_obliquity_ecliptic_IAU(T)
    
    OE_true = 23+26/60+27.407/3600
    
    abs_tol = 1e-6
    
    assert math.isclose(OE, OE_true, abs_tol=abs_tol)
    
    #**************************************************************************
    
    # mean obliquity of the ecliptic (Laskar)
    
    OE = pse.solar.mean_obliquity_ecliptic_Laskar(T)
    
    OE_true = 23+26/60+27.407/3600
    
    abs_tol = 7e-6
    
    assert math.isclose(OE, OE_true, abs_tol=abs_tol)
    
    #**************************************************************************
    
    # nutation in longitude
    
    NL = pse.solar.nutation_longitude(T)
    
    NL_true = -3.788/3600
    
    abs_tol = 2e-6
    
    assert math.isclose(NL, NL_true, abs_tol=abs_tol)
    
    #**************************************************************************
    
    # nutation in obliquity
    
    NO = pse.solar.nutation_obliquity(T)
    
    NO_true = 9.443/3600
    
    abs_tol = 1e-6
    
    assert math.isclose(NO, NO_true, abs_tol=abs_tol)
        
    #**************************************************************************
    
    # true obliquity of the ecliptic
    
    TOE = pse.solar.true_obliquity_ecliptic(T)
    
    TOE_true = 23+26/60+36.850/3600
    
    abs_tol = 1e-6
    
    assert math.isclose(TOE, TOE_true, abs_tol=abs_tol)
    
    #**************************************************************************
    
    # mean anomaly of the sun
    
    MAOTS = pse.solar.mean_longitude_sun(T)
    
    MAOTS_true = -4225.0208 # 94.9792
    
    abs_tol = 1e2
    
    assert math.isclose(MAOTS, MAOTS_true, abs_tol=abs_tol)
    
    #**************************************************************************
    
    # force conditions
    
    error_detected = False    
    try:
        OE = pse.solar.mean_obliquity_ecliptic_Laskar(
            julian_centuries_since_j2000=101)
    except ValueError:
        error_detected = True
    assert error_detected

    #**************************************************************************

#******************************************************************************
#******************************************************************************

def test_julian_day():

    #**************************************************************************
    
    # example 7.a from Meeus (1998), page 61
    
    Y = 1957
    
    M = 10
    
    D = 4.81
    
    calendar = 'gregorian'
    
    JD = pse.solar.julian_day(Y, M, D, calendar=calendar)
    
    true_JD = 2436116.31
    
    abs_tol = 0.01
    
    assert math.isclose(JD, true_JD, abs_tol=abs_tol)
    
    #**************************************************************************
    
    # example 7.b from Meeus (1998), page 61
    
    Y = 333
    
    M = 1
    
    D = 27.5 # +0.5 because it is noon
    
    calendar = 'julian'
    
    JD = pse.solar.julian_day(Y, M, D, calendar=calendar)
    
    true_JD = 1842713.0
    
    abs_tol = 0.1
    
    assert math.isclose(JD, true_JD, abs_tol=abs_tol)
    
    #**************************************************************************
    
    # examples from Meeus (1998), page 62
    
    list_dates = {
        (2000, 1, 1.5, 'gregorian'): 2451545.0,
        (1999, 1, 1.0, 'gregorian'): 2451179.5,
        (1987, 1, 27.0, 'gregorian'): 2446822.5,
        (1987, 6, 19.5, 'gregorian'): 2446966.0,
        (1988, 1, 27.0, 'gregorian'): 2447187.5,
        (1988, 6, 19.5, 'gregorian'): 2447332.0,
        (1900, 1, 1.0, 'gregorian'): 2415020.5,
        (1600, 1, 1.0, 'gregorian'): 2305447.5,
        (1600, 12, 31.0, 'gregorian'): 2305812.5,
        (837, 4, 10.3, 'julian'): 2026871.8,
        (-123, 12, 31.0, 'julian'): 1676496.5,
        (-122, 1, 1.0, 'julian'): 1676497.5,
        (-1000, 7, 12.5, 'julian'): 1356001.0,
        (-1000, 2, 29.0, 'julian'): 1355866.5,
        (-1001, 8, 17.9, 'julian'): 1355671.4,
        (-4712, 1, 1.5, 'julian'): 0.0,
        }
    
    abs_tol = 0.1
    
    for date, true_JD in list_dates.items():
        
        JD = pse.solar.julian_day(date[0], date[1], date[2], calendar=date[3])
                
        assert math.isclose(JD, true_JD, abs_tol=abs_tol)
        
        Y, M, D = pse.solar.calendar_date(JD)
        
        assert math.isclose(Y, date[0], abs_tol=0)
        
        assert math.isclose(M, date[1], abs_tol=0)
        
        assert math.isclose(D, date[2], abs_tol=0)
    
    #**************************************************************************
    
    calendar = 'nocalendar'
    
    try:
        _ = pse.solar.julian_day(
            date[0], 
            date[1], 
            date[2], 
            calendar=calendar)
    except NotImplementedError:
        assert True
    
    #**************************************************************************
    
    error_detected = False    
    try:
        Y, M, D = pse.solar.calendar_date(-1)
    except ValueError:
        error_detected = True
    assert error_detected

    #**************************************************************************
    
    error_detected = False    
    try:
        Y, M, D = pse.solar.calendar_date(-1)
    except ValueError:
        error_detected = True
    assert error_detected
    
    #**************************************************************************
    
#******************************************************************************
#******************************************************************************

def test_calendar_date():
    
    # Source: J. Meeus, Astronomical Algorithms, 1998, page 64, ex. 7.c
    
    # examples from Meeus (1998), page 64
    
    list_dates = {
        2436116.31: (1957, 10, 4.81, 'gregorian'),
        1842713.0: (333, 1, 27.5, 'julian'),
        1507900.13: (-584, 5, 28.63, 'julian'),
        }
    
    abs_tol = 0.1
    
    for JD, true_date in list_dates.items():
    
        Y, M, D = pse.solar.calendar_date(JD)
                
        assert math.isclose(Y, true_date[0], abs_tol=abs_tol)
                
        assert math.isclose(M, true_date[1], abs_tol=abs_tol)
                
        assert math.isclose(D, true_date[2], abs_tol=abs_tol)
    
        JD_new = pse.solar.julian_day(true_date[0], 
                                      true_date[1], 
                                      true_date[2], 
                                      calendar=true_date[3])
                
        assert math.isclose(JD, JD_new, abs_tol=abs_tol)
        
#******************************************************************************
#******************************************************************************

def test_modified_julian_day():

    #**************************************************************************
    
    # example from Meeus (1998), page 63
    
    Y = 1858
    
    M = 11
    
    D = 17
    
    calendar = 'gregorian'
    
    JD = pse.solar.julian_day(Y, M, D, calendar=calendar)
    
    MJD = pse.solar.modified_julian_day(JD)
    
    true_MJD = 0.0
    
    abs_tol = 0.1
            
    assert math.isclose(MJD, true_MJD, abs_tol=abs_tol)
    
#******************************************************************************
#******************************************************************************

def test_julian_day_zero():

    #**************************************************************************
    
    Y = 2000
    
    calendar = 'gregorian'
    
    JD0 = pse.solar.julian_day_zero(Y, calendar=calendar)
    
    # JDE_true = 2451545.0 starts at noon
    
    true_JD0 = 2451543.5 # 
    
    abs_tol = 0.1
            
    assert math.isclose(JD0, true_JD0, abs_tol=abs_tol)
    
    #**************************************************************************
    
    calendar = 'julian'
    
    try:
        JD0 = pse.solar.julian_day_zero(Y, calendar=calendar)
        
    except NotImplementedError:
        assert True
        
    #**************************************************************************

#******************************************************************************
#******************************************************************************

def test_mean_longitude_moon():
    
    JDE = pse.solar.julian_day(1992, 10, 13)
    
    T = pse.solar.time_in_julian_centuries_since_j2000(JDE)
    
    _ = pse.solar.mean_longitude_moon_simple(T)

#******************************************************************************
#******************************************************************************

def test_equation_time_smart1():
    
    # Source: J. Meeus, Astronomical Algorithms, 1998, page 184
    
    # Example 28.a: Find the equation of time on 1992 October 13, at 0h TD
    
    apparent_right_ascension_sun_true = 198.378178 # TODO: calculate this
    nutation_longitude_true = 0.004419
    obliquity_ecliptic_true = 23.4401443
    mean_solar_longitude_true = 201.807193
    
    JDE_true = 2448908.5
    T_true = -0.072183436
    
    true_eq_degrees = 3.427351
    true_eq_minutes = 13.70940
    
    #**************************************************************************
    #**************************************************************************
    
    # spencer
    
    lst = dt.datetime.fromisoformat('1992-10-13 00:00:00+00:00')
    
    doty = pse.timeplace.day_of_the_year(lst)
    
    eq = pse.solar.equation_of_time_minutes_spencer1971(doty)
    
    # minutes ( 15 degrees = 1 hour = 60 minutes )
    
    abs_tol = 0.5
                            
    assert math.isclose(eq, 
                        true_eq_minutes, 
                        abs_tol=abs_tol)
    
    #**************************************************************************
    #**************************************************************************
    
    # time
    
    JDE = pse.solar.julian_day(1992, 10, 13)
    
    abs_tol = 0
                            
    assert math.isclose(JDE, 
                        JDE_true, 
                        abs_tol=abs_tol)
    
    T = pse.solar.time_in_julian_centuries_since_j2000(JDE)
    
    abs_tol = 0
                            
    assert math.isclose(T, 
                        T_true, 
                        abs_tol=abs_tol)
    
    # mean solar longitude
    
    mean_solar_longitude = pse.solar.mean_solar_longitude_vsop87(JDE)
    
    abs_tol = 1e-3
                            
    assert math.isclose(mean_solar_longitude, 
                        mean_solar_longitude_true, 
                        abs_tol=abs_tol)
                        
    # nutation longitude
    
    nutation_longitude = pse.solar.nutation_longitude(T)
    
    abs_tol = 1e-3
                            
    assert math.isclose(nutation_longitude, 
                        nutation_longitude_true, 
                        abs_tol=abs_tol)
    
    # obliquity ecliptic
    
    obliquity_ecliptic = pse.solar.mean_obliquity_ecliptic_IAU(T)
    
    abs_tol = 1e-3
                            
    assert math.isclose(obliquity_ecliptic, 
                        obliquity_ecliptic_true, 
                        abs_tol=abs_tol)
    
    # meeus
    
    eq = pse.solar.equation_of_time_degrees_meeus1998(
        mean_solar_longitude, 
        apparent_right_ascension_sun_true, 
        obliquity_ecliptic, 
        nutation_longitude)
    
    # degrees
    
    abs_tol = 5e-3
                            
    assert math.isclose(eq, true_eq_degrees, abs_tol=abs_tol)

    #**************************************************************************
    #**************************************************************************

#******************************************************************************
#******************************************************************************

def test_equation_time_smart():
    
    # Source: J. Meeus, Astronomical Algorithms, 1998, page 185
    
    # Example 28.b: Find the equation of time on 1992 October 13, at 0h TD
    
    JDE_true = 2448908.5
    T_true = -0.072183436
    earth_orbits_eccentricity_true = 0.016711668
    mean_solar_longitude_true = 201.80720
    obliquity_ecliptic_true = 23.44023 # 23.44023*math.pi/180
    mean_solar_anomaly_true = 278.99397
    
    true_eq_radians = 0.059825572
    true_eq_degrees = 3.427753
    true_eq_minutes = 13+42.6/60
    
    eq =  pse.solar.equation_of_time_radians_Smart(
        mean_solar_anomaly_true*math.pi/180,
        earth_orbits_eccentricity_true,
        mean_solar_longitude_true*math.pi/180,
        obliquity_ecliptic_true*math.pi/180)
    
    # verify
    
    assert math.isclose(eq,
                        true_eq_radians,
                        rel_tol=0.01)
    
    assert math.isclose(eq*180/math.pi,
                        true_eq_degrees,
                        rel_tol=0.01)
    
    # 15 degrees = 1 hour = 60 minutes
                    
    assert math.isclose(eq*(180/math.pi)*60/15,
                        true_eq_minutes,
                        rel_tol=0.01)
    
    #**************************************************************************
    #**************************************************************************
    
    # try to match the inputs, one by one
    
    JDE = pse.solar.julian_day(1992, 10, 13)
                                    
    assert math.isclose(JDE,
                        JDE_true,
                        rel_tol=1e-6)
    
    T = pse.solar.time_in_julian_centuries_since_j2000(JDE)
                                    
    assert math.isclose(T,
                        T_true,
                        rel_tol=1e-6)
                    
    # eccentricity of earth's orbit
    
    eccentricity_earths_orbit = pse.solar.eccentricity_earths_orbit(T)
                                                    
    assert math.isclose(eccentricity_earths_orbit,
                        earth_orbits_eccentricity_true,
                        rel_tol=1e-6)
                    
    # mean obliquity of the ecliptic
    
    OE = pse.solar.mean_obliquity_ecliptic_IAU(T)
                                                    
    assert math.isclose(OE, 
                        obliquity_ecliptic_true,
                        rel_tol=1e-6)
    
    # # true obliquity
    
    # TOE = pse.solar.true_obliquity_ecliptic(T)
                                                                    
    # assert math.isclose(TOE,
    #                     obliquity_ecliptic_true,
    #                     rel_tol=1e-6)
    
    # mean solar longitude
    
    mean_solar_longitude = pse.solar.mean_longitude_sun_simple(T)
                                                                    
    assert math.isclose(mean_solar_longitude,
                        mean_solar_longitude_true,
                        rel_tol=1e-6)
    
    # mean solar anomaly
    
    mean_solar_anomaly = pse.utils.simplify_angle(
        pse.solar.mean_anomaly_sun_simpler(T)
        )+360
    
    # mean_solar_anomaly = pse.solar.mean_anomaly_sun_simpler(T)
                                                                                    
    assert math.isclose(mean_solar_anomaly,
                        mean_solar_anomaly_true,
                        rel_tol=1e-6)
    
#******************************************************************************
#******************************************************************************

# test the incidence angle function

def test_incidence_angle_model():
    
    # Duffie and Beckman (2013): example in page 15
        
    declination_angle = -14*math.pi/180
    latitude = 43*math.pi/180  
    plane_slope_horizontal = 45*math.pi/180
    plane_azimuth_south = 15*math.pi/180
    hour_angle = -22.5*math.pi/180
    
    true_theta = 35*math.pi/180
    
    cos_theta = pse.solar.cosine_incidence_angle(
            declination_angle,  # delta: -23.5 deg <= ... <= 23.5 deg
            latitude,           # phi?
            plane_slope_horizontal,   # beta: -180 <= ... <= 180
            plane_azimuth_south,      # gamma: -180 <= ... <= 180 deg
            hour_angle,         # omega: 
            )
    
    abs_tol = 1*math.pi/180
    
    assert math.isclose(math.acos(cos_theta), true_theta, abs_tol=abs_tol)
    
#******************************************************************************
#******************************************************************************

# test the declination angle function

def test_declination_angle_model():
    
    # Duffie and Beckman (2013): table in page 14
    
    day_list = [
        17,
        47,
        75,
        105,
        135,
        162,
        198,
        228,
        258,
        288,
        318,
        344]
    
    declination_angle_list = [
        -20.9,
        -13.0,
        -2.4,
        9.4,
        18.8,
        23.1,
        21.2,
        13.5,
        2.2,
        -9.6,
        -18.9,
        -23.0]
    
    # this part is not in the example provided
    
    number_sunlight_hours_list = [
        9.216021,
        10.348429,
        11.699134,
        13.186014,
        14.466688,
        15.122768,
        14.824773,
        13.718838,
        12.275835,
        10.79011,
        9.515733,
        8.883039]
    
    latitude = 43
    
    number_decimal_places = 1
        
    # verify each day number and declination angle pair
        
    for (day, true_declination_angle, true_number_sunlight_hours) in zip(
            day_list, declination_angle_list, number_sunlight_hours_list):   
                
        declination_angle = pse.solar.declination_angle(day)
                                                                                                        
        assert math.isclose(round(declination_angle, number_decimal_places), 
                            true_declination_angle)
        
        number_sunlight_hours = pse.solar.number_sunlight_hours(
            declination_angle,
            latitude)
                                                                                        
        assert math.isclose(number_sunlight_hours, 
                            true_number_sunlight_hours,
                            abs_tol=1e-6)
        
        
    #**************************************************************************
    #**************************************************************************
    
#******************************************************************************
#******************************************************************************

# azimuth and zenith angle

def test_azimuth_zenith_angle_model():
    
    # Duffie and Beckman (2013): page 16

    #**************************************************************************
    #**************************************************************************
    
    # test a: latitude 43º, 9:30 AM on February 13
    
    # day_year = 44
    
    latitude = 43*math.pi/180
    
    declination_angle = -14*math.pi/180
    
    hour_angle = -37.5*math.pi/180
    
    true_zenith_angle = 66.5*math.pi/180
    
    true_azimuth_angle = -40.0*math.pi/180 # original: -40.0 (positive west)
    
    # cosine of the zenith angle
    
    abs_tol = 0.1
    
    assert math.isclose(
        pse.solar.cosine_zenith_angle(
            latitude, 
            declination_angle, 
            hour_angle
            ),
        math.cos(true_zenith_angle),
        abs_tol=abs_tol)
    
    # zenith angle
    
    abs_tol = 1*math.pi/180
    
    assert math.isclose(
        math.acos(
            pse.solar.cosine_zenith_angle(
                latitude, 
                declination_angle, 
                hour_angle
                )
            ),
        true_zenith_angle,
        abs_tol=abs_tol)
    
    # azimuth angle
    
    abs_tol = 1*math.pi/180
    
    assert math.isclose(
        pse.solar.solar_azimuth_angle(
            latitude, 
            declination_angle, 
            true_zenith_angle, 
            hour_angle),
        true_azimuth_angle,
        abs_tol=abs_tol)
    
    #**************************************************************************
    #**************************************************************************
    
    # test 2: latitude 43º, 6:30 PM on July 1
    
    # day_year = 182
    
    latitude = 43*math.pi/180
    
    declination_angle = 23.1*math.pi/180
    
    hour_angle = 97.5*math.pi/180
    
    true_zenith_angle = 79.6*math.pi/180
    
    true_azimuth_angle = 112.0*math.pi/180 # original: 112 (positive west)
    
    # cosine of the zenith angle
    
    abs_tol = 0.1
    
    assert math.isclose(
        pse.solar.cosine_zenith_angle(
            latitude, 
            declination_angle, 
            hour_angle
            ),
        math.cos(true_zenith_angle),
        abs_tol=abs_tol)
    
    # zenith angle
    
    abs_tol = 1*math.pi/180
    
    assert math.isclose(
        math.acos(
            pse.solar.cosine_zenith_angle(
                latitude, 
                declination_angle, 
                hour_angle
                )
            ),
        true_zenith_angle,
        abs_tol=abs_tol)
    
    # azimuth angle
    
    abs_tol = 1*math.pi/180
    
    assert math.isclose(
        pse.solar.solar_azimuth_angle(
            latitude, 
            declination_angle, 
            true_zenith_angle, 
            hour_angle),
        true_azimuth_angle,
        abs_tol=abs_tol)
    
    #**************************************************************************
    #**************************************************************************
    
#******************************************************************************
#******************************************************************************

def test_sunrise_sunset_kalogirou():
    
    # Kalogirou, 2013: example 2.6, page 62
    
    # inputs
    
    latitude = 40 # positive north
    #longitude = 26.71875 # negative west (does not matter here)
    time_tolerance = dt.timedelta(minutes=5)
    
    # solution
    
    true_doy = 167
    true_ast_sunrise = dt.datetime.fromisoformat('2013-06-16 04:36:00+02:00')
    true_ast_sunset = dt.datetime.fromisoformat('2013-06-16 19:24:00+02:00')
    
    # calculations
    
    declination_angle = pse.solar.declination_angle(true_doy)

    sunrise_ha = pse.solar.sunrise_hour_angle(declination_angle*math.pi/180, 
                                              latitude*math.pi/180)

    sunset_ha = pse.solar.sunset_hour_angle(declination_angle*math.pi/180, 
                                            latitude*math.pi/180)
    
    ast_sunrise = pse.solar.ast(sunrise_ha, 
                                year=true_ast_sunrise.year, 
                                month=true_ast_sunrise.month,  
                                day=true_ast_sunrise.day, 
                                tzinfo=true_ast_sunrise.tzinfo)

    ast_sunset = pse.solar.ast(sunset_ha, 
                               year=true_ast_sunset.year, 
                               month=true_ast_sunset.month,  
                               day=true_ast_sunset.day, 
                               tzinfo=true_ast_sunset.tzinfo)
    
    # compare
    
    # sunrise
    
    assert ast_sunrise <= true_ast_sunrise + time_tolerance
    
    assert ast_sunrise >= true_ast_sunrise - time_tolerance
    
    # sunset
    
    assert ast_sunset <= true_ast_sunset + time_tolerance
    
    assert ast_sunset >= true_ast_sunset - time_tolerance
    
    #**************************************************************************
    
    # test case without time zone information
    
    latitude = 40 # positive north
    #longitude = 26.71875 # negative west (does not matter here)
    time_tolerance = dt.timedelta(minutes=5)
    
    # solution
    
    true_doy = 167
    true_ast_sunrise = dt.datetime.fromisoformat('2013-06-16 04:36:00')
    true_ast_sunset = dt.datetime.fromisoformat('2013-06-16 19:24:00')
    
    # calculations
    
    declination_angle = pse.solar.declination_angle(true_doy)

    sunrise_ha = pse.solar.sunrise_hour_angle(declination_angle*math.pi/180, 
                                              latitude*math.pi/180)

    sunset_ha = pse.solar.sunset_hour_angle(declination_angle*math.pi/180, 
                                            latitude*math.pi/180)
    
    ast_sunrise = pse.solar.ast(sunrise_ha, 
                                year=true_ast_sunrise.year, 
                                month=true_ast_sunrise.month,  
                                day=true_ast_sunrise.day)

    ast_sunset = pse.solar.ast(sunset_ha, 
                               year=true_ast_sunset.year, 
                               month=true_ast_sunset.month,  
                               day=true_ast_sunset.day)
    
    # compare
    
    # sunrise
    
    assert ast_sunrise <= true_ast_sunrise + time_tolerance
    
    assert ast_sunrise >= true_ast_sunrise - time_tolerance
    
    # sunset
    
    assert ast_sunset <= true_ast_sunset + time_tolerance
    
    assert ast_sunset >= true_ast_sunset - time_tolerance
    
#******************************************************************************
#******************************************************************************

def test_sunrise_sunset_angles():
    
    # Duffie and Beckman (2013): example 1.6.3, page 19
    
    # Calculate the time of sunrise, solar altitude, zenith, solar azimuth, and
    # proﬁle angles for a 60 ◦ sloped surface facing 25 ◦ west of south at 4:00
    # PM solar time on March 16 at a latitude of 43 ◦ . Also calculate the time
    # of sunrise and sunset on the surface.
    
    latitude = 43*math.pi/180
    
    declination_angle = -2.4*math.pi/180
    
    plane_slope_horizontal = 60*math.pi/180
    
    zenith_angle_degrees = 70.3
    
    zenith_angle_radians = round(70.3*math.pi/180,3)
    
    hour_angle = 60*math.pi/180
    
    # results
    
    sunset_hour_angle_degrees = 87.8
    
    sunrise_hour_angle_degrees = -87.8
    
    azimuth_angle_degrees = 66.8
    
    altitude_angle_degrees = 19.7
    
    number_decimal_cases = 1
    
    # assert sunset hour angle
    
    assert (
        round(
            math.degrees(
                pse.solar.sunset_hour_angle(
                    declination_angle, 
                    latitude
                    )
                ),
            number_decimal_cases
            ) == sunset_hour_angle_degrees)
    
    # assert sunrise hour angle
    
    assert (
        round(
            math.degrees(
                pse.solar.sunrise_hour_angle(
                    declination_angle, 
                    latitude
                    )
                ),
            number_decimal_cases
            ) == sunrise_hour_angle_degrees)
    
    # assert solar zenith angle
    
    assert (
        round(
            math.degrees(
                math.acos(
                    pse.solar.cosine_zenith_angle(
                        latitude, 
                        declination_angle, 
                        hour_angle
                        )
                    )
                ),
            number_decimal_cases
            ) == zenith_angle_degrees)

    # assert solar azimuth
    
    assert (
        round(
            math.degrees(
                pse.solar.solar_azimuth_angle(
                    latitude, 
                    declination_angle, 
                    zenith_angle_radians,
                    hour_angle)
                ),
            number_decimal_cases
            ) == azimuth_angle_degrees)   
    
    # assert altitude angle
    
    assert (
        round(
            math.degrees(
                pse.solar.solar_altitude_angle(
                    latitude, 
                    declination_angle, 
                    hour_angle)
                ),
            number_decimal_cases
            ) == altitude_angle_degrees)   
    
    #**************************************************************************
    #**************************************************************************
    
#******************************************************************************
#******************************************************************************

def test_ratios():

    #**************************************************************************
    
    # Duffie and Beckman (2013): example 1.8.1, p. 24, referencing that on p.15  
    
    # Calculate the angle of incidence of beam radiation on a surface located 
    # at Madison, Wisconsin, at 10:30 (solar time) on February 13 if the 
    # surface is tilted 45◦ from the horizontal and pointed 15◦ west of south.
        
    declination_angle = -14*math.pi/180
    latitude = 43*math.pi/180  
    plane_slope_horizontal = 45*math.pi/180
    plane_azimuth_south = 15*math.pi/180
    hour_angle = -22.5*math.pi/180
    
    R_b_true = 1.67
    abs_tol = 0.01
        
    R_b = pse.solar.beam_irradiation_ratio(
        declination_angle, 
        latitude, 
        plane_slope_horizontal, 
        plane_azimuth_south, 
        hour_angle,
        return_average_ratio=False)
        
    #assert round(R_b,number_decimal_cases) == instant_ratio
    
    assert math.isclose(R_b, R_b_true, abs_tol=abs_tol)
    
    #**************************************************************************
    
    # test #2, Example 1.8.2
    
    # Calculate R b for a surface at latitude 40 ◦ N at a tilt 30 ◦ toward the
    # south for the hour 9 to 10 solar time on February 16.
        
    declination_angle = -13*math.pi/180
    latitude = 40*math.pi/180  
    plane_slope_horizontal = 30*math.pi/180
    plane_azimuth_south = 0*math.pi/180
    hour_angle = -45.5*math.pi/180 # 9 to 10 solar time on February 16.
    
    R_b_true = 1.62
    abs_tol = 0.06
        
    R_b = pse.solar.beam_irradiation_ratio(
        declination_angle, 
        latitude, 
        plane_slope_horizontal, 
        plane_azimuth_south, 
        hour_angle,
        return_average_ratio=False)
    
    assert math.isclose(R_b, R_b_true, abs_tol=abs_tol)
    
    #**************************************************************************
    
    # test #3, Example 1.8.3
    
    # Calculate R b for a latitude 40 ◦ N at a tilt of 50 ◦ toward the south 
    # for the hour 9 to 10 solar time on February 16.
        
    declination_angle = -13*math.pi/180
    latitude = 40*math.pi/180  
    plane_slope_horizontal = 50*math.pi/180
    plane_azimuth_south = 0*math.pi/180
    hour_angle = -14*math.pi/180 # 9 to 10 solar time on February 16.
    
    R_b_true = 1.67
    abs_tol = 0.06
        
    R_b = pse.solar.beam_irradiation_ratio(
        declination_angle, 
        latitude, 
        plane_slope_horizontal, 
        plane_azimuth_south, 
        hour_angle,
        return_average_ratio=False)
    
    assert math.isclose(R_b, R_b_true, abs_tol=abs_tol)
    
    #**************************************************************************
    
    # test #4, Example 2.14.1, page 88
    
    # On March 4 at a latitude of 45 ◦ and a surface slope of 60 ◦ determine 
    # R b at 6:30 and R b,ave for the hour 6 to 7 AM .
    
    declination_angle = -7.15*math.pi/180
    latitude = 45*math.pi/180  
    plane_slope_horizontal = 60*math.pi/180
    plane_azimuth_south = 0*math.pi/180
    hour_angle = -82.5*math.pi/180
    
    R_b_true = 44.11 # 39.3
    R_b_avg_true = 4.61 # 4.62
    abs_tol = 0.06
        
    R_b = pse.solar.beam_irradiation_ratio(
        declination_angle, 
        latitude, 
        plane_slope_horizontal, 
        plane_azimuth_south, 
        hour_angle,
        return_average_ratio=False)
    
    assert math.isclose(R_b, R_b_true, abs_tol=abs_tol)
    
    initial_hour_angle = -82.79*math.pi/180
    
    final_hour_angle = -75.0*math.pi/180
    
    R_b_avg = pse.solar.average_beam_irradiation_ratio(
        declination_angle, 
        latitude, 
        plane_slope_horizontal, 
        plane_azimuth_south, 
        initial_hour_angle, 
        final_hour_angle)
    
    assert math.isclose(R_b_avg, R_b_avg_true, abs_tol=abs_tol)
    
    #**************************************************************************
    
#******************************************************************************
#******************************************************************************

def test_extraterrestrial_radiation():
    
    list_examples = [
        # Source: Kalogirou, 2009, page 90, example 2.14
        # March 10th, 2:00 pm solar time at 35°N latitude
        ((35,30),
         'epsg:4326', 
         '2010-03-10 14:00:00+02:00', 
         1383),
        ]
    
    for example in list_examples:
        
        # location
        
        mylocation = pse.timeplace.Location(
            coordinates=example[0], 
            crs_string=example[1])
        
        # ast
        
        myast = dt.datetime.fromisoformat(example[2])
        
        # ast to lst
        
        mylst = pse.solar.ast_to_lst(myast, mylocation)
        
        # day of the year
        
        doty = pse.timeplace.day_of_the_year(mylst)
        
        # extraterrestrial radiation
        
        gsc = pse.solar.extraterrestrial_radiation_normal_duffie2013(doty)
        
        gsc_true = 1383
        
        abs_tol = 1
        
        assert math.isclose(gsc, gsc_true, abs_tol=abs_tol)
        
        gsc = pse.solar.extraterrestrial_radiation_normal_spencer1971(doty)
        
        gsc_true = 1383
        
        abs_tol = 4
        
        assert math.isclose(gsc, gsc_true, abs_tol=abs_tol)
    
#******************************************************************************
#******************************************************************************

def test_air_mass():
    
    # test 1
    
    zenith_angle = 0*math.pi/180
    
    altitude = 0
    
    approx = ''
    
    number_decimal_cases = 1
    
    air_mass = 1
    
    assert math.isclose(
        round(
            pse.solar.get_air_mass(zenith_angle, altitude, approx=approx),
            number_decimal_cases
            ),
        air_mass
        )
    
    #**************************************************************************
        
    # test 2
    
    zenith_angle = 60*math.pi/180
    
    altitude = 0
    
    approx = ''
    
    number_decimal_cases = 1
    
    air_mass = 2
            
    assert math.isclose(
        round(
            pse.solar.get_air_mass(zenith_angle, altitude, approx=approx),
            number_decimal_cases
            ),
        air_mass
        )
    
    #**************************************************************************
        
    # test 3
    
    zenith_angle = 60*math.pi/180
    
    altitude = 0
    
    approx = 'kandy'
    
    number_decimal_cases = 1
    
    air_mass = 2
    
    assert math.isclose(
        round(
            pse.solar.get_air_mass(zenith_angle, altitude, approx=approx),
            number_decimal_cases
            ),
        air_mass
        )
    
#******************************************************************************
#******************************************************************************

# test the lst to ast conversion routine

def test_lst_to_ast():
    
    # tolerance
    
    time_tolerance = dt.timedelta(minutes=1)
    
    list_tests = [
        # format: coordinates, crs, time, dst offset degrees, true standard meridian
        # Example #1: Duffie and Beckman, 2013; ex. 1.5.1, page 11
        # Madison: 43°04′29″N 89°23′03″W
        ((43.0, -89.4),
         'epsg:4326',
         '2010-02-03 10:30:00-06:00', 
         0, 
         -90,
         '2010-02-03 10:19:00-06:00'),
        # Example #2: Kalogirou, 2009; ex. 2.1, page 52
        # Nicosia: 35°10′21″N 33°21′54″E
        ((35.0, 33.3), 
          'epsg:4326', 
          '2010-04-16 10:30:00+02:00', 
          0, 
          30,
          '2010-04-16 10:43:00+02:00'),
        # Example #3: Kalogirou, 2013; ex. 2.3, page 57
        # Athens: 37°59′03″N 23°43′41″E 
        # March 10 at 2:30 pm
        ((37.0+59/60+3/3600, 23+40/60),
          'epsg:4326', 
          '2010-03-10 14:30:00+02:00', 
          0, 
          30,
          '2010-03-10 13:54:00+02:00'), # original: '2010-03-10 14:42:00+02:00'
        ]
    
    for test in list_tests:
        
        mylocation = pse.timeplace.Location(
            coordinates=test[0],
            crs_string=test[1])
        
        mylst = dt.datetime.fromisoformat(test[2])
        
        # compute the local meridian (subtract dst offset since it is included)
        
        local_meridian = mylocation.local_standard_meridian_longitude(
            mylst, test[3])
        
        assert math.isclose(local_meridian,
                            test[4])
        
        # compute the apparent solar time
        
        myast = pse.solar.lst_to_ast(mylst, mylocation, False)
        
        true_ast = dt.datetime.fromisoformat(test[5])
        
        # test upper bound
        
        assert myast <= true_ast + time_tolerance
        
        # test lower bound
        
        assert myast >= true_ast - time_tolerance
        
        # compute the local solar time from the apparent solar time
        
        mylst_from_ast = pse.solar.ast_to_lst(true_ast, mylocation, False)
        
        # test upper bound
        
        assert mylst_from_ast <= mylst + time_tolerance
        
        # test lower bound
        
        assert mylst_from_ast >= mylst - time_tolerance
        
        #**********************************************************************
        
        # compute the solar hour angle using the true ast
        
        solar_ha = pse.solar.hour_angle(true_ast)
        
        mynewast = pse.solar.ast(solar_ha, 
                                 year=true_ast.year, 
                                 month=true_ast.month,
                                 day=true_ast.day,
                                 tzinfo=true_ast.tzinfo)
        
        # compare
        
        # test upper bound
        
        assert mynewast <= true_ast + time_tolerance
        
        # test lower bound
        
        assert mynewast >= true_ast - time_tolerance
        
    #**************************************************************************
    #**************************************************************************
    
    # trigger lst_to_ast and ast_to_lst errors
    
    mylocation.crs_string = 'unknown_crs'
    
    # lst to ast
    
    number_errors = 0
    
    try:
        
        # compute the apparent solar time
        
        myast = pse.solar.lst_to_ast(mylst, mylocation, True)
        
    except ValueError:
        
        number_errors += 1
        
    assert number_errors == 1
    
    # ast to lst
    
    number_errors = 0
    
    try:
        
        # compute the apparent solar time
        
        myast = pse.solar.ast_to_lst(mylst, mylocation, True)
        
    except ValueError:
        
        number_errors += 1
        
    assert number_errors == 1
        
#******************************************************************************
#******************************************************************************